Acknowledgment: This R session is adopted from Ribatet’s short course on “Modelling spatial extremes with the SpatialExtremes package”.

Spatial GEV

Plot the spatial locations

library(SpatialExtremes)
data(wind)
par(mar = rep(0, 4))
maps::map(xlim = c(0, 9), ylim = c(47.5, 56))
points(coord, pch = 15, cex = 0.5)

Fitting spatial GEV

loc.form <- y ~ lon * lat; scale.form <- shape.form <- y ~ 1
#Convert coordinates from degrees to km (rough)
coord[, 1:2] <- 111 * coord[, 1:2]
(M0 <- fitspatgev(wind, scale(coord, scale=FALSE), loc.form, scale.form, shape.form))
##       Model: Spatial GEV model
##    Deviance: 10733.34 
##         TIC: 10766.49 
## 
##     Location Parameters:
##  locCoeff1   locCoeff2   locCoeff3   locCoeff4  
##  2.680e+02  -1.466e-01   1.888e-01   4.548e-05  
##        Scale Parameters:
## scaleCoeff1  
##       33.99  
##        Shape Parameters:
## shapeCoeff1  
##    -0.09664  
## 
## Standard Errors
##   locCoeff1    locCoeff2    locCoeff3    locCoeff4  scaleCoeff1  shapeCoeff1  
##   3.1549723    0.0094921    0.0159158    0.0001288    1.0377212    0.0181140  
## 
## Asymptotic Variance Covariance
##              locCoeff1   locCoeff2   locCoeff3   locCoeff4   scaleCoeff1
## locCoeff1     9.954e+00  -4.245e-03   1.292e-02  -4.676e-05   1.695e+00 
## locCoeff2    -4.245e-03   9.010e-05  -6.018e-05   1.632e-07  -6.348e-04 
## locCoeff3     1.292e-02  -6.018e-05   2.533e-04   7.635e-07   2.467e-03 
## locCoeff4    -4.676e-05   1.632e-07   7.635e-07   1.660e-08  -2.541e-05 
## scaleCoeff1   1.695e+00  -6.348e-04   2.467e-03  -2.541e-05   1.077e+00 
## shapeCoeff1  -3.242e-02   8.809e-05  -6.516e-05   2.652e-07   9.445e-04 
##              shapeCoeff1
## locCoeff1    -3.242e-02 
## locCoeff2     8.809e-05 
## locCoeff3    -6.516e-05 
## locCoeff4     2.652e-07 
## scaleCoeff1   9.445e-04 
## shapeCoeff1   3.281e-04 
## 
## Optimization Information
##   Convergence: successful 
##   Function Evaluations: 547

“Prediction”

library(fields)
x <- seq(min(coord[,1]), max(coord[,1]), length = 100)
y <- seq(min(coord[,2]), max(coord[,2]), length = 100)
grid <- expand.grid(x, y); colnames(grid) <- c("lon", "lat")
## Watch out we scale the data so we need to use the same transformation 
grid[,1] <- grid[, 1] - mean(coord[, 1])
grid[,2] <- grid[, 2] - mean(coord[, 2])
ans <- predict(M0, newdata = grid, ret.per = 20)$Q20
## Do the plot
maps::map(xlim = range(x) / 111, ylim = range(y) / 111)
image(x / 111, y / 111, matrix(ans, 100), add = TRUE, col = tim.colors(64))
contour(x / 111, y / 111, matrix(ans, 100), add = TRUE); maps::map(add = TRUE)

Bayesian Hierarchical Model

Plot the spatial domain and station locations

data(rainfall)
data(swissalt)
par(mar = rep(0, 4), ps = 16)
image(lon.vec, lat.vec, alt.mat, col = terrain.colors(64), asp = 1,
      bty = "n", xlab = "n", ylab = "n", axes = FALSE)
swiss(add = TRUE, city = TRUE)
points(coord, pch = 16, cex = 0.75, col = "blue")

Prior specification

hyper <- list()
hyper$betaMeans <- list(loc = rep(0, 3), scale = rep(0, 3), shape = 0)
hyper$betaIcov <- list(loc = diag(rep(1/1000, 3)), scale = diag(rep(1 / 1000, 3)), shape = 1 / 10)
hyper$sills <- list(loc = c(1, 12), scale = c(1, 1), shape = c(1, 0.04))
hyper$ranges <- list(loc = c(5, 3), scale = c(5, 3), shape = c(5, 3))
hyper$smooths <- list(loc = c(1, 1), scale = c(1, 1), shape = c(1, 1))
prop <- list(gev = c(3, 0.1, 0.3), ranges = c(1, 0.8, 1.2), smooths = rep(0, 3))
start <- list(sills = c(10, 10, 0.5), ranges = c(20, 10, 10), smooths = c(1, 1, 1), beta = list(loc = c(25, 0, 0), scale = c(33, 0, 0), shape = 0.001))

Running the Gibbs sampler

loc.form <- scale.form <- y ~ lon + lat; shape.form <- y ~ 1
chain <- latent(rain, coord[, 1:2], "powexp", loc.form, scale.form, shape.form,
                hyper = hyper, prop = prop,
start = start, n = 100, burn.in = 500, thin = 5)
chain
## Effective length: 100 
##          Burn-in: 500 
##         Thinning: 5 
##    Effective NoP: 75.52989 
##              DIC: 29166.1 
## 
##   Regression Parameters:
##       Location Parameters:
##            lm1       lm2       lm3     
## ci.lower    1.09608  -0.01467  -0.19594
## post.mean  30.50024   0.03854  -0.12421
## ci.upper   64.82665   0.08884  -0.05800
## 
##          Scale Parameters:
##            lm1       lm2       lm3     
## ci.lower   -5.30524   0.00193  -0.05280
## post.mean   5.95500   0.01644  -0.03197
## ci.upper   17.89916   0.03166  -0.01415
## 
##          Shape Parameters:
##            lm1    
## ci.lower   0.06993
## post.mean  0.16673
## ci.upper   0.26259
## 
## 
##   Latent Parameters:
##       Location Parameters:
##         Covariance family: powexp 
##            sill    range   smooth
## ci.lower    5.718  13.520   1.000
## post.mean   9.930  23.864   1.000
## ci.upper   17.044  45.431   1.000
## 
##       Scale Parameters:
##      Covariance family: powexp 
##            sill     range    smooth 
## ci.lower    0.1274   5.9886   1.0000
## post.mean   0.3122  17.6629   1.0000
## ci.upper    0.5957  27.9217   1.0000
## 
##       Shape Parameters:
##      Covariance family: powexp 
##            sill       range      smooth   
## ci.lower    0.004151  13.704008   1.000000
## post.mean   0.007827  23.628722   1.000000
## ci.upper    0.013391  36.571825   1.000000

Making 20-year return level map

map.latent(chain, ret.per = 20, plot.contour = FALSE, col = tim.colors())

Max-stable processes

USHCN data

data(USHCNTemp)
maps::map("usa")
points(metadata[, c("lon","lat")], pch = 16, cex = 0.75)

F– madogram

data <- maxima.summer
coord <- as.matrix(metadata[, c("lon", "lat")])
fmadogram(data, coord)

fmadogram(data, coord, which = 'ext', n.bins = 100, pch = 16)
abline(h = 1, lty = 2); abline(h = 2, lty = 2)

Simulation from various max-stable models

x <- seq(0, 10, length = 100)
n.sim <- 5
schlather <- rmaxstab(n.sim, x, "powexp", nugget = 0, range = 3, smooth = 1)
extremalt <- rmaxstab(n.sim, x, "tpowexp", DoF = 4, nugget = 0, range = 3, smooth = 1)
brown <- rmaxstab(n.sim, x, "brown", range = 3, smooth = 1)
smith <- rmaxstab(n.sim, x, "gauss", var = 2)

plot(x, schlather[1,], type = "l", ylim = range(schlather))
for (i in 2:5) lines(x, schlather[i,], col = i)

plot(x, smith[1,], type = "l", ylim = range(smith))
for (i in 2:5) lines(x, smith[i,], col = i)

Two-step max-stable process fitting

data(rainfall)
alt <- coord[,3]; coord <- coord[,-3]
## Transformation to unit Frechet margins
rain.frech <- apply(rain, 2, gev2frech, emp = TRUE)
(fit <- fitmaxstab(rain.frech, coord, "powexp", nugget = 0))
##         Estimator: MPLE 
##             Model: Schlather 
##          Weighted: FALSE 
##    Pair. Deviance: 1136875 
##               TIC: 1137468 
## Covariance Family: Powered Exponential 
## 
## Estimates
##   Marginal Parameters:
##   Assuming unit Frechet.
## 
##   Dependence Parameters:
##   range   smooth  
## 38.4402   0.8528  
## 
## Standard Errors
##  range  smooth  
##  8.713   0.119  
## 
## Asymptotic Variance Covariance
##         range     smooth  
## range   75.90902  -0.91383
## smooth  -0.91383   0.01417
## 
## Optimization Information
##   Convergence: successful 
##   Function Evaluations: 71

One-step model fitting

loc.form <- y ~ lon + lat + alt
scale.form <- y ~ lon + lat
shape.form <- y ~ 1
(fit <- fitmaxstab(rain, coord, "powexp", nugget=0, loc.form, scale.form, shape.form,
marg.cov = cbind(alt = alt)))
## Computing appropriate starting values
## Starting values are defined
## Starting values are:
##        range       smooth    locCoeff1    locCoeff2    locCoeff3    locCoeff4 
## 27.257344723  0.974059614 29.764826111  0.037985247 -0.135302350  0.008372902 
##  scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1 
## 10.529521199  0.013058818 -0.040111724  0.156381541
##         Estimator: MPLE 
##             Model: Schlather 
##          Weighted: FALSE 
##    Pair. Deviance: 2232232 
##               TIC: 2243533 
## Covariance Family: Powered Exponential 
## 
## Estimates
##   Marginal Parameters:
##     Location Parameters:
## locCoeff1  locCoeff2  locCoeff3  locCoeff4  
## 28.660445   0.039237  -0.134847   0.008425  
##        Scale Parameters:
## scaleCoeff1  scaleCoeff2  scaleCoeff3  
##     8.59677      0.01565     -0.03937  
##        Shape Parameters:
## shapeCoeff1  
##      0.1911  
##   Dependence Parameters:
##   range   smooth  
## 29.7692   0.9751  
## 
## Standard Errors
##       range       smooth    locCoeff1    locCoeff2    locCoeff3    locCoeff4  
##   9.6080167    0.1744583    6.9550328    0.0097462    0.0166063    0.0005552  
## scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1  
##   4.8427858    0.0071586    0.0115508    0.0518712  
## 
## Asymptotic Variance Covariance
##              range       smooth      locCoeff1   locCoeff2   locCoeff3 
## range         9.231e+01  -1.482e+00   6.902e+00   1.389e-03  -1.886e-02
## smooth       -1.482e+00   3.044e-02  -6.333e-02   1.543e-05   9.781e-05
## locCoeff1     6.902e+00  -6.333e-02   4.837e+01  -5.307e-02  -4.118e-02
## locCoeff2     1.389e-03   1.543e-05  -5.307e-02   9.499e-05  -4.612e-05
## locCoeff3    -1.886e-02   9.781e-05  -4.118e-02  -4.612e-05   2.758e-04
## locCoeff4     7.735e-04  -1.761e-05   5.940e-04  -1.975e-06   2.533e-06
## scaleCoeff1   3.595e+00  -5.025e-02   1.111e+01  -1.342e-02  -3.458e-03
## scaleCoeff2   9.077e-03  -5.956e-05  -1.033e-02   2.306e-05  -2.246e-05
## scaleCoeff3  -1.761e-02   8.477e-05  -1.183e-02  -1.005e-05   6.792e-05
## shapeCoeff1   2.473e-01  -2.315e-03   4.184e-02  -2.615e-05  -9.515e-05
##              locCoeff4   scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1
## range         7.735e-04   3.595e+00    9.077e-03   -1.761e-02    2.473e-01 
## smooth       -1.761e-05  -5.025e-02   -5.956e-05    8.477e-05   -2.315e-03 
## locCoeff1     5.940e-04   1.111e+01   -1.033e-02   -1.183e-02    4.184e-02 
## locCoeff2    -1.975e-06  -1.342e-02    2.306e-05   -1.005e-05   -2.615e-05 
## locCoeff3     2.533e-06  -3.458e-03   -2.246e-05    6.792e-05   -9.515e-05 
## locCoeff4     3.082e-07  -4.400e-04    1.178e-06   -1.332e-06   -6.617e-06 
## scaleCoeff1  -4.400e-04   2.345e+01   -2.758e-02   -1.355e-02    4.509e-02 
## scaleCoeff2   1.178e-06  -2.758e-02    5.125e-05   -3.135e-05    2.846e-07 
## scaleCoeff3  -1.332e-06  -1.355e-02   -3.135e-05    1.334e-04   -1.169e-04 
## shapeCoeff1  -6.617e-06   4.509e-02    2.846e-07   -1.169e-04    2.691e-03 
## 
## Optimization Information
##   Convergence: successful 
##   Function Evaluations: 911
plot(fit)

LRT Model selection

M0 <- fitmaxstab(rain, coord, "powexp", nugget = 0, locCoeff4 = 0, loc.form, scale.form,
                 shape.form, marg.cov = cbind(alt = alt))
## Computing appropriate starting values
## Starting values are defined
## Starting values are:
##        range       smooth    locCoeff1    locCoeff2    locCoeff3    locCoeff4 
## 27.257344723  0.974059614 29.764826111  0.037985247 -0.135302350  0.008372902 
##  scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1 
## 10.529521199  0.013058818 -0.040111724  0.156381541
M1 <- fitmaxstab(rain, coord, "powexp", nugget = 0,
loc.form, scale.form, shape.form, marg.cov = cbind(alt = alt))
## Computing appropriate starting values
## Starting values are defined
## Starting values are:
##        range       smooth    locCoeff1    locCoeff2    locCoeff3    locCoeff4 
## 27.257344723  0.974059614 29.764826111  0.037985247 -0.135302350  0.008372902 
##  scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1 
## 10.529521199  0.013058818 -0.040111724  0.156381541
anova(M0, M1)
## Eigenvalue(s): 161.12 
## 
## Analysis of Variance Table
##    MDf Deviance Df Chisq Pr(> sum lambda Chisq)    
## M0   9  2246053                                    
## M1  10  2232232  1 13821              < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model selection via TIC

M1 <- fitmaxstab(rain, coord, "powexp", nugget = 0, loc.form, scale.form, shape.form,
                 marg.cov = cbind(alt = alt))
## Computing appropriate starting values
## Starting values are defined
## Starting values are:
##        range       smooth    locCoeff1    locCoeff2    locCoeff3    locCoeff4 
## 27.257344723  0.974059614 29.764826111  0.037985247 -0.135302350  0.008372902 
##  scaleCoeff1  scaleCoeff2  scaleCoeff3  shapeCoeff1 
## 10.529521199  0.013058818 -0.040111724  0.156381541
M2 <- fitmaxstab(rain, coord, "whitmat", loc.form, scale.form,
                 shape.form, marg.cov = cbind(alt = alt))
## Computing appropriate starting values
## Starting values are defined
## Starting values are:
##        nugget         range        smooth     locCoeff1     locCoeff2 
##  9.754459e-05  3.511644e+01  3.877117e-01  2.976483e+01  3.798525e-02 
##     locCoeff3     locCoeff4   scaleCoeff1   scaleCoeff2   scaleCoeff3 
## -1.353023e-01  8.372902e-03  1.052952e+01  1.305882e-02 -4.011172e-02 
##   shapeCoeff1 
##  1.563815e-01
TIC(M1, M2)
##      M1      M2 
## 2243533 2243665

Parameters and return level maps

xlim <- range(M2$coord[,1]); ylim <- range(M2$coord[,2])
data(swissalt)##Retrieve altitude at other places
idx.x <- which(lon.vec >= xlim[1] & lon.vec <= xlim[2])
idx.y <- which(lat.vec >= ylim[1] & lat.vec <= ylim[2])
x <- lon.vec[idx.x]; y <- lat.vec[idx.y]
n.x <- length(x); n.y <- length(y)
covariates <- array(alt.mat[idx.x, idx.y], dim = c(n.x, n.y, 1), dimnames = list(NULL, NULL, "alt"))
covariates[is.na(covariates)] <- 0
par(mar = rep(0, 4), mfrow = c(1, 3))
SpatialExtremes::map(M2, x, y, covariates, param = "loc",
                     xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)
SpatialExtremes::map(M2, x, y, covariates, param = "scale",
                     xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)
SpatialExtremes::map(M2, x, y, covariates, param = "quant", ret.per = 20,
                     xaxt = "n", yaxt = "n", col = tim.colors()); swiss(add = TRUE)